Influence of polydispersity on micromechanics of granular materials 
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We study the effect of polydispersity on tiie macroscopic physical properties of granular packings 
in two and three dimensions. A mean-field approach is developed to approximate the macroscale 
quantities as functions of the microscopic ones. We show that the trace of the fabric and stress 
tensors are proportional to the mean packing properties (e.g. packing fraction, average coordination 
number, and average normal force) and dimensionless correction factors, which depend only on 
the moments of the particle-size distribution. Similar results are obtained for the elements of the 
stiffness tensor of isotropic packings in the linear afline response regime. Our theoretical predictions 
are in good agreement with the simulation results. 
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I. INTRODUCTION 



The physics of granular media has received a lot of at- 
tention because of its scientific challenges and industrial 
relevance. The structural and dynamical properties of 
granular materials differ from those of ordinary solids, 
liquids, or gases due to nonlinearity and disorder [ll-Q- 
On the microscopic level, a static assembly of grains con- 
sists of particles which interact with their neighbors in 
order to prevent interpenetration. In spite of the uni- 
form density of granular packings, the resulting contact 
and force networks between particles are highly inhomo- 
geneous (^-lal, leading to many intriguing phenomena in 
these systems. Describing the behavior via micromechan- 
ical approaches, in which the discrete nature of the sys- 
tem is taken into account, is thus commonly preferred 
to continuum-mechanical approaches where some heuris- 
tic assumptions have to be made in order to construct 
the constitutive equations for macroscopic fields. One 
can then express the macroscopic physical quantities in 
terms of the microscale ones. For example, thermal and 
electrical conductivities are related to the trace of the 
fabric tensor, a micro-geometrical probability of the ori- 
entations of contacts. While the relationship between 
macroscopic and microscopic properties of granular me- 
dia has been studied widely |l|,l3|,|7[, the question remains 
as to what extent the macroscale quantities are sensitive 
to the micro-scale details, and how large is the error in- 
troduced in the calculation of the "observable quantities" 
by taking into account only the average packing proper- 
ties. 
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Granular materials in nature and industry consist of 
particles with the common property of polydispersity. 
It is known that size polydispersity affects the mechan- 
ical behavior of granular systems (e.g. shear strength) 
P, 0] as well as their space-filling properties (e.g. pack- 
ing fraction) [l3, lll[ , which are crucial in many chemical 
processes like absorption, filtering, etc. Polydispersity 
in most studies, so far, has been restricted to narrow 
size distributions mainly to prevent long-range structural 
order; however, there are a few studies where broader 
ranges of particle-size distribution are investigated pj|.lll|- 
[l3|. In this paper, we address the question of how devi- 
ation from the monodisperse case influences the macro- 
scopic properties of granular assemblies. 

We consider a special case of spherical particles [or 
disks in two dimensions (2D)] allowing for analytical cal- 
culations. The main goal is to develop a mean-field ap- 
proach to calculate the desired microscopic quantities, 
such as the trace of the fabric and stress tensors, and 
the elements of the stiffness tensor in two- and three- 
dimensional polydisperse granular systems. These quan- 
tities are directly connected to macroscopic quantities 
such as thermal and electrical conductivities, isotropic 
pressure, and bulk and shear moduli. A similar ana- 
lytical approach has been already used in Ref. [ij] to 
calculate the trace of the fabric tensor in 2D packings, 
where it turned out that the trace of fabric is factor- 
ized into three contributions: (i) the volume fraction, (ii) 
the mean coordination number, and (iii) a dimensionless 
correction factor which only depends on the particle-size 
distribution. Using a similar approach, here we investi- 
gate also the stress and stiffness tensors and extend the 
method to 3D cases. In order to compare the analyti- 
cal results with numerical simulations, we first construct 
static packings of grains using contact dynamics simula- 
tions |15l - [l7| . The initial dilute systems of rigid parti- 
cles are compressed by imposing a confining pressure to 



get the final static homogeneous packings [18]. Compar- 
isons are then made between the results of our mean-field 
model and the exact values obtained from the numerical 
simulations. 

This work is organized in the following manner: The 
fabric tensor of a polydisperse assembly of spherical par- 
ticles is investigated in Sec.|TTl and a mean-field approach 
is introduced to calculate the trace of fabric. We present 
the analytical results for the calculation of the stress ten- 
sor in Sec. IIII[ and the same approach is used in Sec. |IV] 
to investigate the stiffness tensor elements in frictionless 
isotropic packings. In Sec. |Vl the analytical calculations 
are compared to numerical simulations of corresponding 
packings of polydisperse particles. Finally, we discuss 
and conclude the results in Sec. lVIl Detailed calculations 
for two-dimensional packings of disks are presented in the 
Appendix. 




FIG. 1: Schematic picture showing a typical particle with 
radius a surrounded by identical particles of average radius 
(a) in a 3D packing of spheres. 



II. FABRIC TENSOR 

A. Single-particle case 

Various definitions of the fabric tensor have been used 
in the literature to describe the spatial arrangement of 
the particles in a granular assembly |19l - [2l| . The fabric 
tensor of the second order for one particle is defined as 
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where Cp is the number of contacts of particle p, and l„ 
is the a component of the branch vector I , connect- 
ing the center of particle p to its contact c. In the case 
of spherical particles, the unit branch vector l"" /\l'"'\ 
and the unit normal vector n at contact c are identi- 
cal. The trace of the single-particle fabric tensor in a 
i'-dimensional system is 
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i.e. the number of contacts of particle p. 
B. Many- particle case 



The average fabric tensor (/i^^) ^ enables us to describe 
the global contact network in a given volume V . Assum- 
ing that the contribution of particle p (lying inside V) 
to the average fabric tensor is proportional to its volume 
Vp, we obtain 



1 "^ 



(3) 



where the sum runs over all particles lying inside V ^ 
and (• • ■)y denotes the volume weighted average. Us- 
ing Eq. (0) to calculate the trace of the average fabric 
tensor, we get 
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which can be interpreted as the contact number den- 
sity. Alternative possibilities, e.g. using the volume of 
the polygon that contains the particle (obtained e.g. via 
Voronoi tessellation), or introducing constant prefactors 
or slightly different volume contributions are not dis- 
cussed here (see Refs. [H, HO, [13, HI for more details). 
In a monodisperse packing, Eq. dH for identical particles 
is reduced to {h^^)^ =<^^, where is the packing fraction 
(<j) — ^ Vp/V), and z is the average coordination num- 
ber {z — '^pCp/N). We note that only "real" contacts 
contribute to the calculation of z, and geometrical neigh- 
bors without a permanent physical contact, which do not 
contribute in the fabric and force carrying structures, are 
not considered here. 



C. Polydispersity 

For an accurate evaluation of the trace of the aver- 
age fabric tensor in a polydisperse granular packing, one 
should take into account the contributions from all par- 
ticles. However, if the distribution function of particle 
radii is known, {h^^)^ can be approximated as a func- 
tion of the moments of the size distribution. We assume 
a polydisperse distribution of particle radii with proba- 
bility f{a)da to find the radius between a and a + da, and 
with L f(a)da = 1. The continuum limit of Eq. Q is 
then given by 
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Here, C{a) is the average coordination number of parti- 
cles with radius a. We evaluate C{a) usin g a mean-field 
approach similar to the one proposed in [2^ and used 
already in [ij] to study the trace of the fabric tensor. 
In the following, we concentrate on the case of spher- 
ical particles in three-dimensional systems (see the de- 
tailed calculations for two-dimensional packings of disks 
in the Appendix). Let us suppose that each particle in 
the polydisperse granular medium is surrounded by iden- 
tical particles of average radius (a) (see Fig. [T]), where 

(a) = / af{a)da. The surface of a reference particle of 
radius a is then shielded by its C{a) neighboring particles 
of radius (a) . The space angle covered by a neighboring 
particle on the reference particle in a three-dimensional 
packing of spheres is 



n{a) =27r 1 
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(6) 



The total fraction of shielded surface, also called linear 
compacity, is obtained as 

C{a) 

cJa) = ^ V n{a)a^ = Vt{a)Cia)lATT. (7) 

Now, another basic assumption is that the total fraction 
of shielded surface Cs is independent of the particle radius 
a. As a result, the expected mean coordination number 
becomes 



C{a)f{a)da = iTTCsq^, 



(8) 



with go=/o I{a)l^{a)da. Using Eqs. © and ^ one 
finds 



C{a) = 



q„n{a)' 



(9) 



The trace of the fabric tensor for a polydisperse packing 
is then obtained by substitution of Eq. © in Eq. ([5]), 



where the correction factor g^ is defined as 
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Here, {a^) and (a''') denote the fcth moments of the 
size distribution /(a) and the modified distribution 
f{a)/V,{a) normalized by q^, respectively. We note that 
5 J depends only on the size distribution function /(a). 



D. Narrow size distributions 

By introducing e(a)=a/(a) — 1, which ranges between 
-1 and oo depending on the choice of a, Eq. ^ can be 




1/a (exact) 

1/il (first order approximation) 

1/ii (second order approximation) 
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FIG. 2: l/n(a) as a function of e. The exact value (solid line) 
is compared with the first-order (dashed line) and second- 
order (dash-dotted line) approximations. The inset shows 
more clearly the deviation of the approximations from the 
exact value. 



written as 



n(a) 



27r 1 
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Indeed, e(a) quantifies the deviation from the mean par- 
ticle size (a) [e.g. e(a) equals zero in the monodisperse 
case]. Hence, for narrow size distributions, we approx- 
imate l/ri(a) by Taylor expansion around e—0 (corre- 
sponding to Taylor expansion around a={a)). By Taylor 
expansion to the second order in e, one obtains 
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and Cj 



(2-V3)7r' 1 2x73(2-73)2 

The first-order approximation deviates 



3(3+y3)(2-y3)27r 

significantly from the exact value (see Fig. [2]). However, 
the second-order expansion provides a good approxima- 
tion with less than 1% error in the range —0.5 < e < 7.5 
(or 0.5(a) <a<8. 5(a)). 

Therefore, the correction factor [Eq. ([TT]) ] for narrow 
size distributions becomes 
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Equation P^ should account for arbitrarily shaped size 
distributions /(a) as long as they are not too wide. Note 
the different nomenclature in Ref. [3J|, where the above 
equation is introduced with different abbreviations and 
coefficients. 
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FIG. 3: (a) A typical contact c between the reference particle 
p and its neighboring particle, (b) The contact unit vectors 
n*'^ if, and if^ (c) The normal [F^") and tangential [Ff") 
components of the contact force F . 



III. STRESS TENSOR 



For a spherical grain, we project the contact unit vectors 
(n , i , t ) onto an arbitrary Cartesian coordinate sys- 
tem [Figs, mja) andlS^b)], and write the force-averaged 
stress tensor of a single particle as 
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where the Wmnki function is defined as 

Wmnki = sin"(6'c) cos"(6'c) sin''(.^c) cos\ipc), 



(17) 



(18) 



with 0^9c<TT and 0^ipc<2TT. Using Eq. (|T6|). the trace 
of the stress tensor becomes 
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Equation (IT^ remains valid also in the 2D case (see Ap- 
pendix). As expected for isotropic packings, the trace of 
the stress tensor and therefore the isotropic pressure P 
(=o'q,q/3) do not depend on the tangential forces. 



A. Single-particle case 

The micromechanical expressions for the components 
of the stress tensor cr^^ of a single particle in a static 
granular assembly are [23, [27| 
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where F is the force exerted on particle p by its neigh- 
boring particle at contact c. 

One could assume in a crude approximation that the 
force at contact c is equal to PPfiP'^ + Ff t^'^ + Ff t^'^ in 
a three-dimensional system, where F^, F[^, and F^^ are 
the average normal and tangential contact forces around 
the particle p, and n*"^, if, and if are the normal and 
tangential unit vectors at contact c, respectively. Then 
the force- averaged stress tensor becomes 
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B. Many-particle case 

In the many-particle case, the average stress tensor in 
a given volume V is defined as [23| 

N N Cp 

(0.-^E^.-:.-^EEC>r. (20) 
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where the sum runs over all particles lying inside V . Us- 
ing Eq. p9)) to calculate the trace of the average stress 
tensor, we get 

N N 

i^..)v = -i7 E ^P^'o^c. = T; E ""P^nCp. (21) 
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C. Polydispersity 

Now we assume a polydisperse distribution of particle 
radii with probability f(a)da to find the radius between 

a and a + da, and with J f{a)da = 1. Assuming that 



the average contact force exerted on a particle depends 
only on its radius a, the continuous limit of Eq. (j2ip in 
a mean-field approximation is given by 






aFn{a)C{a)f{a)da. 



(22) 



In Eq. (im, it is supposed that all particles of size a 
have a certain mean coordination number C{a) and a 
certain mean normal force Fn{a). Indeed, particles of 
the same size may have different coordination number 
and normal contact forces, however, the main goal here is 
to propose a method to calculate macroscopic quantities 
without taking into account all the microscopic details of 
the system. We use the mean-field approach introduced 
in Sec. Ill CI to evaluate C{a). By substitution of Eq. © 
in Eq. (I22|) we get 

,^ , N ir aFn{a){^da 
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According to the mean-field approach used in Sec. Ill CI 
C{a) increases with increasing radius a. Now, let us as- 
sume that the average normal force Fn{a) also increases 
with a, so that the ratio Fn{a)/C{a) remains roughly con- 
stant [28]. We calculate this ratio for the average-sized 
particles in the following: 
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therefore 
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By substitution of Eq. ([25)) in Eq. ([23|) , we obtain 
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D. Narrow size distributions 
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In the limit of narrow size distributions, we approxi- 
mate 1/Vi^{a) by Taylor expansion around e = (similar 
toSecUD]): 
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with A^= A 
1 



B^ 



and C, 
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4(2-v^)4.2 - 18(2-^)3.2 - By substitution of Eq. 

Eq. (P?]) we obtain the correction factor g^ for arbitrary 

narrow distributions, 
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IV. STIFFNESS TENSOR 
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The linear response of a material to "weak" external 
perturbations is described by a fourth rank tensor, which 
is called the elastic or stiffness tensor [3, [23] . This ten- 
sor has 81 and 16 elements in three- and two-dimensional 
systems, respectively, but they are not all independent. 
Symmetry considerations reduce the number of indepen- 
dent elements. For example, the elastic behavior of 
isotropic materials can be described by only two inde- 
pendent parameters, usually represented by Lame coef- 
ficients A and ii. In this section, the stiffness tensor of 
a homogeneous and isotropic assembly of polydisperse 
particles is investigated (for the case of an anisotropic 
monodisperse system, see, e.g. J30l.l3l|). 

The stiffness tensor for a spherical particle, where 
affine deformation is assumed, is defined as [3, [S^j 

2a2 '^''' 

^a,p,-i,ri~~r7' 2^^ "-^a '^^ ^■y "?) + ktU^ tp U^ t^ ), 
^ c— 1 

(30) 
where t'^'^ is the unit vector parallel to the tangential 
component of the contact force F'P'^ [see Fig. [3l^c)]. The 
volume weighted average of C is then given by 

1 ^ 

N Cp 

V E 2alY.{KK^nfnP^nP^+hK%%?^^tP/). (31) 

p—l c—l 

Note that the stiffness tensor is basically determined by 
the packing geometry. For ease of calculation, we con- 
sider only frictionless packings, i.e. kt is set to zero here- 
after. Using the microscopic information of the contact 
orientations, one can accurately calculate the elements of 
C via Eq. ((3T|) . Next, the Lame constants /i and A can 
be deduced from the stiffness tensor, e.g. as ^—{C^^^^)^ 
and A-|-2/i=(Cjjjj)^,, or, more generally, as A=(C....)j, and 



A+2^=(C,^,J^, where 
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and D is the dimension of the system. The macroscopic 
physical quantities of interest are the bulk modulus K 
and the shear modulus G, which can be deduced from 
the Lame coefficients in isotropic materials as 



and for narrow size distributions, one obtains 
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and 



^,, ,,^2 (q„.)v+(^-l)(q.„). ,„,. 
X/fc„=(A+— /.t)/fc„== — '-^ — . (34) 

Now, assuming a polydisperse probability distribution 
of particle radii /(a), Eq. (PT|) for kt=0 can be written 
as 
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Since the packings are supposed to be isotropic and ho- 
mogeneous, we assume that grains are scattered homo- 
geneously around the reference particle. Therefore, the 
summation over neighbors can be approximated by the 
following integration in three dimensions (for the 2D case, 
see Appendix): 
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We present the reduced form of the fourth rank tensor 
by mapping 0/3(777) -^ i{j), i.e. 11 -> 1, 22 -^- 2, 33 —> 3, 
12 ^ 4, 13 ^ 5, and 23 -^ 6. Using Eqs. ®, ([33, and 
((36)) , one obtains 
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where W^, elements were defined in Eq. p^ . After 
integration on and iy9, the volume weighted average of 
the stiffness tensor for an isotropic polydisperse packing 
becomes 
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The correction factor g,j^ is defined as 
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with the same coefficients as defined after Eqs. p^ and 
(1281) . To summarize this section, the Lame constants for 
frictionless isotropic packings are 



fi = \^ {kn(t)zg.^) I (107r(a)), 
and the shear and bulk moduli are 

G/fc„ = (^z53)/(107r(a)), 
and 

ir/fc„ = (0z53)/(67r(a)). 



(41) 
(42) 
(43) 



Notably, KjG = 5/3 in three-dimensional frictionless 
isotropic packings, independent of their size distribution 
and average packing properties. 



V. SIMULATION RESULTS 

To verify the theoretical predictions of the previous 
sections, we carry out numerical simulations with the 
help of the contact dynamics (CD) algorithm [ISl - llTj . We 
first construct 2D and 3D static homogeneous packings in 
zero gravity by compressing the initial dilute configura- 
tion of particles [Fig. |4] (left)]. Periodic boundary condi- 
tions are imposed in all directions to avoid the side effects 
of lateral walls. The compaction is achieved by impos- 
ing a constant external pressure Pext and letting the size 
of the system evolve in time [33] ■ As the volume of the 
system decreases, after a while particles touch each other 
and build an inner pressure Pinn, which resists and even- 
tually compensates Pext, so that finally Pinn equals Pext- 
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FIG. 4: Schematic of a 2D granular system subjected to a 
constant external pressure: the initial dilute gas (left), and 
the final homogeneous packing (right). Periodic boundaries 
are marked with dashed lines. 



TABLE I: Properties of three different types of polydisperse 
packings generated with uniform size distributions, w denotes 
the width of each distribution (w = o,max~^min)- 



type symbol a^^,, 

SMPl • 0.67 

SMP2 D 0.40 

SMP3 A 0.22 



"max 'Jmax/"„ 

1.34 2 

1.60 4 

1.76 8 



w/2{a) {a')/{af 
0.34 1.04 

0.60 1.12 

0.77 1.19 



Particles prevent further compaction, and a static ho- 
mogeneous configuration is reached [Fig. |4] (right)]. The 
fuU description of the packing generation method can be 
found in [l8|. In order to ihustrate the vahdity range of 
our assumptions, we generate three types of polydisperse 
packings with uniform particle-size distributions but with 
different widths (see Table |l|. We denote the samples 
SMPl, SMP2, and SMP3, respectively, by fuU circles, 
open squares, and full triangles throughout this section. 
To investigate the effect of friction, we construct a new 
packing for each value of the particle-particle friction co- 
efficient fi, . In particular, the results corresponding to 
l^f=0, 0.1, and 1.0 are hereafter denoted by green, blue 
and red colors. The number of grains contained by pack- 
ings are 3000 and 10000 in 2D and 3D cases, respectively. 
For comparison with the theory, we first test the va- 
lidity of assumptions made in Sec. Ill CI The linear com- 
pacity Cs is displayed in Fig. \S\ for the static configura- 
tions of particles obtained from the isotropic compression 
simulations. For each particle p, the surface angle fl^ 
covered by its neighboring particle at contact c is calcu- 
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FIG. 5: (Color online) Linear compacity Cs as a function 
of particle radius a for two-dimensional (top) and three- 
dimensional (bottom) packings constructed with different size 
distribution widths and different friction coefficients p,. 
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FIG. 6: (Color online) Contact number C(a) as a function 
of particle radius a for two-dimensional packings. The lines 
correspond to the mean-field approximation of C{a) according 
to Eq. (O. 



lated, and the linear compacity of particle p is obtained 
as Cg = ^^Z;^ ^c /^TT or Cg = ^^Z;^ rij, /47r for two- or three- 
dimensional packings, respectively. Next, we divide the 
range of possible values of the particle radius a into 25 
bins. Each data point in Fig. [5] corresponds to the mean 
value of Cs, averaged over all particles in the same bin. 
The contribution of the rattler particles, which transmit 
no force, is excluded. For moderate widths of size distri- 
butions (SMPl), Cs is approximately constant in a for a 
given packing (we note that the fiuctuations of Cg around 
its mean value in a given packing originate from the finite 
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FIG. 7: (Color online) The same plots as in Fig. [6l but for 
three-dimensional packings. 
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FIG. 8: The estimated value of the trace of the average fabric tensor (j>zg-^ vs the exact value {h^^)^ obtained from the 
simulations, for several (a) two- and (b) three-dimensional samples. The dashed lines indicate the identity. The insets show 
more precisely that the deviations increase with w, but remain less than 5% in all cases. The symbols are chosen the same as 
in Table [H 



size of the samples). However, Cs is remarkably above 
the average vahie for small particle sizes in wider dis- 
tributions (SMP2 and SMP3). This is a common prop- 
erty of our highly polydisperse packings (with uniform 
size distribution) that the fraction of shielded surface is 
larger than the average for small particles if rattlers are 
excluded (see [ij| for uniform volume distributions). A 
similar behavior has been observed in discrete element 
method simulations of soft particles [34|. There, it is 
also shown that if rattlers are included in the statistics, 
the small particles on average are less covered than the 
larger ones. However, the deviation of small particles 
from the average Cs decreases as the volume fraction of 
the packing increases by incremental compression. 

Another point is that Cs depends strongly on the di- 
mension of the system and the friction coefficient. In- 
creasing the friction fi, stabilizes the system in a less 
dense state and decreases the connectivity of the contact 
network [35|, |36| . Therefore, we expect lower values of Cg 
and C{a) when increasing /i^ , as confirmed by the data. 

In Figs. [6] and [3 the coordination number C{a) is 
shown as a function of a for the same set of systems 
as in Fig. [5l For comparison, we also plot C{a) from 
Eq. ^ . Here, the average coordination number z of the 
packing is taken from the simulation results, n{a) is pro- 
vided by Eq. ^ or (|Aip , and the size distribution of each 
packing after the compaction process is used to calculate 
Qg. The mean-field approach of Sec. Ill CI qualitative Iv 
fits well to the data, however, the slopes of the curves 
are slightly greater than the corresponding slopes of the 
best- fit curves over the data points (not shown). Con- 
sequently, one expects that the mean-field approach to 
calculate the trace of the fabric tensor {h^J)^ leads to 
somewhat overestimated values. For each packing, we 
calculate the exact value of {h^^)y via Eq. ^ and com- 
pare it with the mean-field approximation [Eq. (1101) ]. Fig- 



ure [5] reveals that Eq. (|10l) slightly overestimates {h^J)^ 
in both two- and three-dimensional systems. The devi- 
ation increases with the width of the size distribution, 
but remains less than 5% in all cases. For comparison, 
note that g^ can reach up to 1.19 and 1.45 in 2D and 
3D uniform samples, respectively (see Table HI)) : there- 
fore, ignoring the correction factor would cause up to 
19% and 45% error, respectively. 

Next, we investigate the average properties of the con- 
tact force network. In Sec. IIII CI we applied the mean- 
field approach of Sec. Ill CI to estimate the isotropic pres- 
sure in a given polydisperse granular sample. However, 
due to the presence of the normal component of the con- 
tact force Fn{a) in Eq. 1^^, one needs to make one fur- 
ther assumption about the particle-size dependence of 
Fn{a) to be able to calculate the integral and obtain 
(paa)v from the average quantities. 

The simulation results [Fig. EKa)] reveal that the aver- 
age normal force exerted on the particle is an increasing 
function of the particle radius for 2D and 3D (the con- 
tribution of rattlers is again excluded) . With increasing 
friction coefficient and w, the average normal force in- 



TABLE II: Correction factors in two and three dimensions 
for uniform size distributions SMPl, SMP2, and SMP3 intro- 
duced in Table in 
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FIG. 9: (Color online) (a) Normal component of the contact force Fn{a), averaged over all particle radii in the same bin, in 
terms of the particle radius a for two-dimensional (top) and three-dimensional (bottom) packings constructed with different 
size distribution widths and different friction coefficients fi.. (b) F„{a) scaled by the contact number C(o) for the same set of 
samples as in (a). 



creases. This is reminiscent of the behavior of C{a) as a 
function of a (Figs.[n]and[71). Interestingly, the increasing 
rates are similar in both figures. Therefore, it is reason- 
able to assume that the ratio F„(a)/Cia) is independent 
of a, as already observed in 2D [2^. Figure[SlJb) confirms 
the validity of this assumption. We note that the fiuctua- 
tions in Fig.[9][b) are reduced as the system size increases. 
In Fig.[Tni we compare the exact value of {(Jaa)v '^itli the 
corresponding value from Eq. ([25]) [or Eq. (|A10p ]. which 
is obtained based on the above assumption. The results 
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FIG. 10: The estimated value of the trace of the stress ten- 
sor using Eq. (|A10|l (top), and Eq. ([26|l (bottom) divided by 
{o'aa)v obtained directly from the simulations. Each data 
point corresponds to a different 2D (top) or 3D (bottom) 
packing using symbols as in Table ID 



are in reasonable agreement with theory for both two- 
and three-dimensional packings, with a standard devia- 
tion of 2% to 6% for increasing w. 



(a) 



0.31 



CM 



0.30 
K 0.29 
0.28 
0.27 
0.26 



0.15 




/- 






□ ,-''' 


0.14 




_,/' 




•/' 


- 


0.13 




- 



0.26 



0.28 

K/k„ 



0.30 



0.13 



0.14 
G/k„ 



0.15 



(t") 0.100 






0.062 


• / 




0.095 

K 

CD 


/ 

u / 


A 

o 


0.058 


a / 




'"co 0.090 
0.085 






0.054 
0.050 








0.085 0.090 0.095 0.100 




0.050 0.054 0.058 0.06 


2 




K/k„ 






G/k„ 





FIG. 11: (Color online) The estimated values of the bulk K 
(left) and shear G (right) moduli according to Eqs. (|42p and 
(|43)) in 3D [Eqs. (|A18|) and (1A19P in 2D] vs the values obtained 
from the simulation results. Each data point corresponds to 
one frictionless sample and the dashed lines indicate the iden- 
tity. The results are separately shown for (a) 2D and (b) 3D 
samples. The same symbols as in Table |T] are used. 
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Finally, we turn to the calculation of the stiffness ten- 
sor elements for isotropic materials. We note that to 
evaluate the true elastic moduli, one should apply an 
incremental strain and measure the resulting change of 
the stress tensor. Alternatively, one can read the moduli 
from the elements of the stiffness tensor, assuming the 
affine motion of the particles, which cannot be taken for 
granted however, and which is the subject of future stud- 
ies. Here, using the packing configuration obtained from 
the simulation, we calculate the elements of the average 
stiffness tensor via Eq. pip . Next, the elastic moduli of 
the packing are calculated using Eqs. (15^ . ([55]) and ([M|. 
The results are then compared to the estimated values of 
the bulk and shear moduli calculated via Eqs. (|42l) and 
(|43l) [or Eqs. (|A18|) and (|A19P ]. Figure [n] displays the 
results for several two- and three-dimensional packings; 
the agreement is satisfactory within a 5% error (also in 
the case of frictional packings which is not shown here) . 

According to our analytical results, the ratio between 
the bulk and shear moduli K/G is 5/3 for isotropic pack- 
ings independent of z, 0, and even the size distribution. 
This suggests that in isotropic packings, the ratio be- 
tween the P-wave velocity Vp—J{K+^G)/p and the 

5- wave velocity Vg^^jG/ p is always \/3. An experimen- 
tal test shows that Vp/Vg for a compressed polydisperse 
packing of glass beads remains around 1.7 over a wide 
range of pressures from 1 to 7 MPa [33| (see also |38|'). 
Note, however, that anisotropic regular lattice structures 
do not necessarily show the same ratio |31| . 



VI. DISCUSSION AND CONCLUSION 

In conclusion, a mean-field approach is developed to 
isolate the influence of size polydispersity on the phys- 
ical properties of granular assemblies. We are inter- 
ested in how the microscale quantities are linked to the 
macroscale ones. 

We find that the trace of fabric and stress tensors fac- 
torize into the mean packing properties (for example, av- 
erage coordination number, packing fraction, and aver- 
age normal contact force) and dimensionless correction 
factors, which depend on the moments of the particle- 
size distribution (and approach unity for monodisperse 
packings). The method is extended to estimate the ele- 
ments Q ,,, of the stiffness tensor. This tensor describes 
the linear afRne response of the packing to weak exter- 
nal perturbations, when practically the contact network 
between the particles remains unchanged. The elements 
C ,,, are also proportional to the average quantities and 
a dimensionless correction factor, which is a function of 
the size distribution. 

Numerical simulations illustrate the validity range of 
our analytical predictions and of the assumptions on 
which the mean-field method is based. We note that the 
deviation of the macroscopic quantities of interest from 
the average packing properties increases with increasing 
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FIG. 12: The dimensionless correction factors g^ in terms 
of the width w of the uniform size distribution in (a) two 
and (b) three dimensions. ■w/2{a) = corresponds to the 
monodisperse case. 



the width w of the particle-size distribution. Figure [12] 
shows the summarized correction factors g- as a function 
of the width u" of a uniform size distribution, with the 
average particle size (a). Neglecting the correction fac- 
tors would cause remarkable errors, especially for wide 
distributions. Interestingly, g^ is insensitive to the width 
of the size distribution in the 3D case. Therefore, ac- 
cording to Eqs. P^ and ((33]), we expect that the elas- 
tic moduli of a polydisperse packing of spheres is only 
moderately affected by the choice of w. The results of 
molecular dynamics (MD) simulations of soft frictionless 
spheres imply (see Eq. (12) in [3^) that the bulk modu- 
lus does not depend on the width of the size distribution, 
in agreement with our analytical results. 

The predictive value of this mean-field method should 
be examined also by comparing the theoretical predic- 
tions with experimental data. For a direct comparison, 
one needs to measure the average packing properties, 
e.g. z and (f), which are not easily accessible in exper- 
iments (even though microcomputed tomography (Mi- 
croCT) scan determines the geometry with micrometer 
accuracy nowadays [39|). Alternatively, by elimination 
of 0z between our analytical results, one obtains linear 
relationships between the macroscopic physical proper- 
ties via some coefficients, which depend on the moments 
of the size distribution. Such linear relations between 
macroscopic quantities have been investigated in the lit- 
erature, e.g. between the elastic moduli and conductivity 
[40| or isotropic pressure [41j, and can be verified ex- 
perimentally. Future studies will more closely examine 
the nonaffinity of deformations of isotropic as well as 
anisotropic packings of frictional and possibly even co- 
hesive particles. 
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Appendix A: Analytical results in two dimensions 

Fabric tensor . In a two-dimensional packing of disks 
[Fig. [131a)], the surface angle covered by a neighboring 
particle on the reference particle is 



r2(a) — 2arcsin 



(Al) 



a+ (a) I' 
and the total fraction of the shielded surface is given by 

^ C{a) 

c,{a) = - — V n(a)a = n(a)C(a)/27r. (A2) 

4=1 

Assuming that Cs is independent of a, one can write the 
mean coordination number z as 



(A3) 



C{a)f{a)da = 2TrCsqo- 



Equations (U2|) and (IM]) lead again to Eqs. ^ and (HH 
for C{a) and (/i„„)vi with the correction factor 



V{a)i}^da , 

il[a) _ {a ) 



(A4) 



Qo / V{a)f{a)da 
Jo 

By introducing e=a/{a) — l, we rewrite Eq. (JA1|) as 

1 \ 



fl{a) =2 arcsin 



2 + e ' 



(A5) 



and approximate l/f2(a) to the first order in e for narrow 
size distributions 



1 



n{a) 



A[+B[e, 



(A6) 



where A[ = ^ and B[ = ^^ . Figure \W[h) reveals that 
the approximation has a less than 1% error in the range 
-0.5 < e < 1.3 (or 0.5(a) < a < 2.3(a)). Hence, g^ for 
narrow size distributions becomes 



B' 



5i 



K \{a){a') 



(A7) 




Ma (exact) 

1/i2 (first order approximation) 
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FIG. 13: (a) A typical particle with radius a surrounded by 
identical particles of average radius (a) in a 2D packing of 
disks. The thick solid arcs show the shielded surface of the 
central particle, (b) l/il(a) as a function of e in two dimen- 



Stress tensor . For a two-dimensional disk, by disre- 
garding the z direction, i.e., in the x—y plane [by requiring 
61 = f and Ff = {) in Fig.lHb)], one obtains 



5^=^ 
K 



ppy^i cos2(v?) sin((^) cos((^)| 



+FPy^i~^'^^^^^'^°^^^^ cos2((^) 
*'^V -sin2('/') sin(</3)cos(v?)^ 



(A8) 



and its trace 



Cp D 



^P c=l a=l ^ 

c 

Using Eqs. (P5)) and (P5|) . the average stress tensor in 2D 
becomes 



i^.J. 



TT{a) 



(AlO) 



with 



52 



Jo ^ (g) 



da 



(All) 



By Taylor expansion around e 
l/f72(a) as 



0, we approximate 



1 



172(a) 



with A' = A' = ^ and B' = 

2 1 j^2 2 

approximated by 



A' +B'e 

2 2 

18\/3 



(A12) 



Therefore, g^ can be 



52 



2_ 

A' 



1- 



KM 



A'J{a^y 



(A13) 



Stiffness tensor . Similarly to the three-dimensional 
analysis presented in Sec. IIVI we approximate 
the summation over neighbors in Eq. (|35p by 

—^ J nanpn^riridO, which leads to the following 
reduced stiffness tensor (by mapping 11 -^- 1, 22 ~> 2 
and 12 ^- 3): 



(C), 




(A14) 



with 



,(=g,) = {a')J{a^), 
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(A15) 



which for narrow size distributions is approximated as 



B' 






-1 



(A16) 



In two dimensions, one finds that the Lame constants for 
frictionless isotropic packings are 



fj, = X= {kn(t)zg^) / (An), 



and, hence, the shear and bulk moduli are 



and 



G/kn = {^zg,)/{A7r), 



K/kn = {(bzg^)/{2TT). 



(A17) 



(A18) 



(A19) 
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